10. 2要因実験の分析
https://gyazo.com/4f7dfd91756d579cff6ea820ee6ff329
10.1. 2つの要因とは
「知覚時間」の実験を、少しやり方を変えてS1氏が測定した
前章と同じように要因$ Aは測定の「条件」
$ a=4であることも20回ずつ測定していることも変わらない
ただし、「条件」の4つの水準は2つの異なった「状態」に区分されている
「平常時」に10回、「起床直後」に10回
1要因の実験デザインは実験計画法の基礎であり、頻繁に利用される重要なモデル
しかし測定値の高低に影響する要因は、1つ取り上げれば、それでいつでも十分というわけにはいかない
その場合は2つの目の要因を導入する
今回は要因$ Bとして「状態」を導入する
水準数は$ b = 2
要因Aと要因Bの水準の組み合わせによって表現される区分をセル(cell)という ここには要因Aの4水準×要因Bの2水準で8つのセルがある
セル$ jk: 要因$ Aの$ j番目の水準と要因$ Bの$ k番目の水準
table: 表10-2 S1氏のデータの平均と標準偏差
平常時 起床直後
対照 聴音 音読 運動 対照 聴音 音読 運動
平均 30.08 32.21 32.07 29.46 29.54 31.83 29.48 31.80
標準偏差 2.27 2.43 1.52 2.04 1.79 2.70 2.61 2.10
https://gyazo.com/3fd8ebc7d531191511f4700eaa60bbaf
10.1.1. モデル構成
$ y_{ijk} = \mu + a_j + b_k + (ab)_{jk} + e, \quad e \sim N(0, \sigma_e) \qquad (10.1)
左辺の$ y_{ijk}はセル$ jkにおける$ i番目の測定値
例えば$ y_{321} = 30.64
要因の影響を受ける前の測定値の平均
右辺第2項$ a_jは要因$ Aの水準$ jの主効果
たとえば、要因$ A「条件」にもし主効果があるならば、要因$ B「状態」の水準によらず、「対照」「聴音」「音読」「運動」のどれかの間に一定の差があるということ
右辺第3項$ b_kは要因$ Bの水準$ kの主効果
たとえば、要因$ B「状態」にもし主効果があるならば、要因$ A「条件」の水準によらず、「平常時」と「起床直後」の測定値に一定の差があるということ
右辺第4項$ (ab)_{jk}は要因$ ABのセル$ jkの交互作用効果
交互作用をあらわすカッコは掛け算でないことを表現している
一方の要因の水準の違いによって、他方の要因の水準の間の平均のパタンが異なる状態
主効果だけでは説明がつかない場合、「要因Aと要因Bには交互作用がある」という
右辺第5項の$ eはセル内の散らばりを表現しており、誤差変数
$ eはセルによらず平均$ 0、標準偏差$ \sigma_eの正規分布に従うことが仮定される
10.1.2. モデルの制約
(9.10)式に準じて要因Aの水準の効果の和は$ 0である
ただし$ a_aだけ残して$ a_1から$ a_{a-1}は右辺に移行し、ここでは
$ a_a = (-1) \times (a_1 + \cdots + a_{a-1}) \qquad (10.2)
なぜこのように変形するかというと、右辺の$ a-1個が自由に推定できる母数であり、左辺の$ a_aは生成量であることを明示するため
要因Bの効果の和も$ 0であり、同様の表記で
$ b_b = (-1) \times (b_1 + \cdots + b_{b-1}) \qquad (10.3)
要因Bのなかで自由に推定できる母数の数は$ b-1個であり、左辺の$ b_bは生成量とする
交互作用は添え字$ iと$ jの両方に関して和が$ 0であるという制約を入れ
$ \begin{aligned} (ab)_{ak} & = (-1) \times ((ab)_{1k} + (ab)_{2k} + \cdots (ab)_{(a-1)k}), \\ (ab)_{jb} & = (-1) \times ((ab)_{j1} + (ab)_{j2} + \cdots (ab)_{j(b-1)}) \qquad (10.4) \end{aligned}
したがって、交互作用ABのなかで自由に推定できる母数の数は$ (a-1)(b-1)個
$ \muと$ \sigma_eの2個も母数であるから、独立した2要因計画のモデルには、合計で$ ab + 1個($ =(a-1)+(b-1)+(a-1)(b-1)+2)の自由に推定できる母数があることがわかる
表10-1でいうならば、(10.1)式中の母数は$ 9個($ =4\times2+1)
$ \bm \theta = (a_1, a_2, a_3, b_1, (ab)_{11}, (ab)_{21}, (ab)_{31}, \mu, \sigma_e) \qquad (10.5)
(10.1)式中の生成量は
$ (a_4, b_2, (ab)_{12}, (ab)_{22}, (ab)_{32}, (ab)_{42}, (ab)_{41}) \qquad (10.6)
モデル式(10.1)式を書き直す
$ \mu_{jk} = \mu + a_j + b_k + (ab)_{jk} \qquad (10.7)
$ y_{ijk} = \mu_{jk} + e, \quad e \sim N(0, \sigma_e) \qquad (10.8)
(10.7)式はセル$ jkの平均を表現している
(10.8)式を観察すると、セルの平均にセル内の誤差が加わって測定値が生成されていることがわかる
(10.8)の左式右辺第1項に$ ab個の平均があり、右式に$ 1個の標準偏差があることから、自由に推定できる母数の数は$ ab+1個であることを再確認できる
10.1.3. 事後分布
(10.8)式より、測定値の確率分布は以下の正規分布である
$ f(y_{ijk}|\mu_{jk}, \sigma_e) \qquad (10.9)
セル$ jk内で$ n_{jk}個の測定が独立されているならば、セル内の測定値$ \bm y_{jk} = (y_{1jk}, \cdots, y_{ijk}, \cdots, y_{n_{jk}jk})の同時確率分布は
$ f(\bm y_{jk}|\mu_{jk}, \sigma_e) = f(y_{1jk}|\mu_{jk}, \sigma_e) \times \cdots f(y_{ijk}|\mu_{jk}, \sigma_e) \times f(y_{n_{jk}jk}|\mu_{jk}, \sigma_e) \qquad (10.10)
セルごとに測定値の数値は異なっていて構わない
データ全体を$ \bm y = (\bm y_{11}, \cdots \bm y_{1b}, \bm y_{21}, \cdots\cdots, \bm y_{ab})と表記し、セルごとの平均をまとめて$ \mu = (\mu_{11}, \cdots, \mu_{1b}, \mu_{21}, \cdots\cdots, \mu_{ab})と表記すると(2.12)式に相当する尤度は以下になる
$ f(\bm y|\bm\theta) = f(\bm y|\bm\mu, \sigma_e) = f(\bm y_{11}|\mu_{11}, \sigma_e) \times \cdots \times f(\bm y_{ab}|\mu_{ab}, \sigma_e) \qquad (10.11)
ここで母数ベクトル$ \bm\thetaは$ (\bm\mu, \sigma_e)でもよいのであるが、以下のように表記したほうが、あとの分析過程が明快になる
$ \bm\theta = (\mu, \bm a, \bm b, (\bm{ab}), \sigma_e) \qquad (10.12)
ただし、$ \bm a = (a_1, \cdots, a_{a-1}), \bm b = (b_1, \cdots, b_{b-1}), (\bm{ab}) = ((ab)_{11}, \cdots, (ab)_{1\ b-1}, (ab)_{21}, \cdots\cdots, (ab)_{a-1\ b-1})である
$ \bm\thetaのなかの$ \muと$ \sigma_eの事前分布として、ここでは次のように仮定した
$ \mu \sim U(0, 100), \quad \sigma_e \sim U(0, 50) \qquad (10.13)
主効果$ \bm a, \bm bと交互作用効果($ \bm{ab})に関しては特に範囲を定めず、無限と言って構わない程の広い一様分布とした
(2.14)式に相当する同時事前分布を
$ \begin{aligned} f(\bm\theta) & = f(\mu) \times f(a_1) \times \cdots \times f(a_{a-1}) \times f(b_1) \times \cdots \times f(b_{b-1}) \\ & \times f((ab)_{11}) \times \cdots \times f((ab)_{1\ b-1}) \times \cdots \times f((ab)_{a-1\ b-1}) \times f(\sigma_e) \end{aligned}
とし、(2.15)式に相当する事後分布を導く
$ f(\bm\theta|\bm y) \propto f(\bm y|\bm\theta)f(\bm\theta) \qquad (10.14)
MCMC法により、事後分布・生成量・予測分布に従う乱数を生成することが可能
10.2. 推測例1(交互作用の分析)
$ 21000個の乱数を5本発生させ、バーンイン期間を$ 1000とし、$ T=100000の乱数によって母数の事後分布を近似した
MCMCの設定に関しては、後述する推測例2, 3でも共通
10.2.1. 母数の推定値
母数の推定結果を表10-3に示す
https://gyazo.com/a76e9e9e0e1eb3796a52e8d7a1095020
$ b_2 = -b_1であり$ (ab)_{j2} = -(ab)_{j1}なので省略している
主効果や交互作用は、それに含まれる項の少なとも1つが$ 0でないときに効果があると判定する
95%の確信で判定するならば、$ a_1は高々$ -0.24であるし、$ a_2は少なくとも$ 0.45なので、母数の推定結果から要因Aの効果はあると判定できる
また$ (ab)_{31}は少なくとも$ 0.39であり、$ (ab)_{41}は高々$ -0.55なので、交互作用効果もあると判定できる
10.2.2. 水準とセルの効果の評価
主効果があるとしたらどの水準が(ドメイン知識に照らして)基準$ cより大きい(あるいは小さい)のだろうか
交互作用があるとしたらどのセルが基準$ cより大きいのだろうか
これらに関するもう1つの判定法を紹介する
「研究仮説:水準・交互作用の効果は$ cより大きい」が正しい確率は、生成量(9.12)式と以下の生成量のEAPで評価する
$ u_{b_k>c}^{(t)} \begin{cases} 1 & b_k^{(t)} > c \\ 0 & それ以外の場合 \end{cases} \quad u_{(ab)_{jk}>c}^{(t)} = \begin{cases} 1 & (ab)_{jk}^{(t)} > c \\ 0 & それ以外の場合 \end{cases}
ここでは入門的分析ということで$ c = 0として計算した確率を表10-4に示す
table: 表10-4 水準・交互作用の効果が0より大きい(小さい)確率 (S1)
a₁ a₂ a₃ a₄ b₁ (ab)₁₁ (ab)₂₁ (ab)₃₁ (ab)₄₁
0 < 0.016 0.995 0.472 0.347 0.712 0.604 0.538 0.993 0.003
≦ 0 0.984 0.005 0.528 0.653 0.288 0.396 0.462 0.007 0.997
$ a_1は$ 98.4\%で$ 0以下
$ a_2は$ 99.5\%で$ 0より大きい
$ (ab)_{31}は$ 99.3\%で$ 0より大きい
$ (ab)_{41}は$ 99.7\%で$ 0以下
以上から要因Aの趣向かと交互作用があると判定できる
10.2.3. 要因の効果の評価
主効果と交互作用効果の有無ではなく、その効果の大きさはどれ程だろうか
(10.1)式右辺の第1項以外の4つの項が互いに独立であるとすると測定値の分散は単純な和となる
$ \sigma_y^2 = \sigma_a^2 + \sigma_b^2 + \sigma_{ab}^2 + \sigma_e^2 \qquad (10.15)
ここで
$ \begin{aligned} \sigma_a^2 & = \frac{1}{a}(a_1^2 + \cdots + a_a^2), \quad \sigma_b^2 = \frac{1}{b}(b_1^2 + \cdots + b_b^2), \\ \sigma_{ab}^2 & = \frac{1}{a \times b}((a \times b)個の(ab)_{jk}^2の総和) \end{aligned}
これは確率変数の性質であるから、水準ごとのデータ数が異なっても影響されない
要因の効果の大きさを解釈するために利用できる1つの指標としては、説明率がある
$ \eta_a^2 = \frac{\sigma_a^2}{\sigma_y^2}, \quad \eta_b^2 = \frac{\sigma_b^2}{\sigma_y^2}, \quad \eta_{ab}^2 = \frac{\sigma_{ab}^2}{\sigma_y^2}, \quad \eta_t^2 = \frac{\sigma_a^2 + \sigma_b^2 + \sigma_{ab}^2}{\sigma_y^2}
説明率は測定値の分散に占める、要因の分散の比だった
$ \eta_t^2は2つの要因と交互作用の分散の和による説明率
要因の効果の大きさを解釈するために利用できるもう一つの指標としては、効果量がある
$ \delta_a = \frac{\sigma_a}{\sigma_e}, \quad \delta_b = \frac{\sigma_b}{\sigma_e}, \quad \delta_{ab} = \frac{\sigma_{ab}}{\sigma_e} \qquad (10.16)
効果の大きさに関する生成量の推定結果を表10-5に示す
https://gyazo.com/712151d9af1b348e4607978fe1ac38b5
交互作用ABの効果の標準偏差は$ 0.956(0.254)[0.469, 1.464] であり、およそ$ 0.96秒
交互作用ABの説明率は$ 0.127(0.057)[0.031, 0.250] であり、約1割3分
交互作用ABの効果料は$ 0.405(0.110)[0.193, 0.624] であり、$ 40.5\%
10.2.4. セル平均の事後分布
交互作用があると判定された場合には、主効果の有無によらず交互作用の分析を進める
S1氏のデータには、交互作用と要因Aの主効果の両方があったので、交互作用の分析を進める
一方の要因の水準の違いによって、他方の要因の水準の間の平均の高低のパタンが異なるならば、主効果の解釈がしにくいため
2要因の分析で交互作用効果の存在が確信されたら、一方の要因の水準ごとの他方の要因の水準間の差を推測するとデータに対する理解が深まる
まずは(10.7)式を利用してセル$ jk平均の生成量を計算し、事後分布を近似する
$ \mu_{jk}^{(t)} = \mu^{(t)} + a_j^{(t)} + b_k^{(t)} + (ab)_{jk}^{(t)} \qquad (10.17)
結果を表10-6に示す
https://gyazo.com/17345f884157abd19b821b19c1819de1
EAP推定値は表10-2の標本平均と実質的に同じ
「研究仮説:セル平均$ \mu_{jk}は$ \mu_{j'k'}より大きい」が正しい確率は、以下のEAPで評価する
$ u_{\mu_{jk}>\mu_{j'k'}}^{(t)} = \begin{cases} 1 & \mu_{jk}^{(t)} > \mu_{j'k'}^{(t)} \\ 0 & それ以外の場合 \end{cases} \qquad (10.18)
10.2.5. 特に興味のある2セル間の比較
ここでは平均の差、効果量、非重複度、優越率、閾上率を選んで例示する
標準偏差には前節までに求めた$ \sigma_e^{(t)}を利用する
要因「条件」の水準「音読」における、他方の要因「状態」の「平常時」と「起床直後」間の差の推測を行い、結果を表10-7に示した
https://gyazo.com/eacab8219363eba12083ecdf633f7f95
閾上率の差には$ 1.0秒を指定した
要因「状態」の水準「起床直後」における、他方の要因「条件」の「運動」と「音読」間の差の推測を行い、結果を表10-8に示した
https://gyazo.com/f2aafe9801c9398ba8411bf55b2c8754
10.3. 推測例2(3水準以上の主効果の分析)
本節と次節では主効果のみがある場合の分析例を示す
前節とは別人のS2氏が自身で測定した「知覚時間」のデータが表10-9
図10-2に平均値プロットを示し、表10-10に水準・交互作用の効果が$ 0より大きい、または小さい確率を示した
https://gyazo.com/828b5a21e5144e1cd3e82849965cd925
table: 表10-10 水準・交互作用の効果が0より大きい(小さい)確率 (S2)
a₁ a₂ a₃ a₄ b₁ (ab)₁₁ (ab)₂₁ (ab)₃₁ (ab)₄₁
0 < 0.207 0.992 0.994 0.000 0.358 0.781 0.490 0.496 0.228
≦ 0 0.793 0.008 0.006 1.000 0.642 0.219 0.510 0.504 0.772
$ a_2が$ 0より大きい確率は$ 0.992
$ a_3が$ 0より大きい確率は$ 0.994
$ a_4が$ 0以下である確率は有効数字3桁で$ 1.000
以上のことから、要因Aの主効果だけがあると判定される
平均値プロットを観察すると、2本の折れ線がよく似た形状で、同じ高低を示している
これが要因Aの主効果だけがある場合の特徴
S2氏のデータだけを分析する時は表10-3や表10-5に相当する表も作成する
しかしここでは表10-5のなかの要因Aに相当する効果の大きさに関する指標を示す
table: 表10-11 効果の大きさに関する生成量の推定結果 (S2)
EAP post.sd 2.5% 5% 50% 95% 97.5%
ηₐ² 0.227 0.071 0.092 0.111 0.225 0.347 0.369
δₐ 0.558 0.117 0.330 0.366 0.557 0.753 0.790
説明率は$ 0.227(0.071)[0.092, 0.369]
要因Aの主効果はデータの変動の$ 22.7\%を説明している
効果量は$ 0.558(0.117)[0.330, 0.790]
条件が一定に保たれているセル内での散らばりの$ 55.8\%の変動を要因Aの主効果が生み出している
10.3.1. 水準間の比較
(10.18)式にならって「研究仮説:水準の効果$ a_jは$ a_{j'}より大きい」が正しい確率を以下のEAPで評価する
$ u_{\mu_{a_j}>\mu_{a_{j'}}}^{(t)} = \begin{cases} 1 & a_{j}^{(t)} > a_{j'}^{(t)} \\ 0 & それ以外の場合 \end{cases}
table: 表10-12 行jの水準の効果が列j'の水準の効果より大きい確率
条件 対照(a_1) 聴音(a_2) 音読(a_3) 運動(a_4)
対照(a_1) 0.000 0.024 0.021 0.979
聴音(a_2) 0.976 0.000 0.476 1.000
音読(a_3) 0.979 0.524 0.000 1.000
運動(a_4) 0.021 0.000 0.000 0.000
「聴音」は「対照」より$ 0.976の確率で大きい
「音読」は「対照」より$ 0.979の確率で大きい
「対照」は「運動」より$ 0.979の確率で大きい
10.3.2. 連言命題が正しい確率
表10-12の考察より、要因Aの水準は『(「聴音」「音読」)>「対照」>「運動」』の3つの群に分かれそう
連言命題で表現される研究仮説「$ (a_2, a_3) > a_1 > a_4が正しい確率を考える
この研究仮説が正しい確率は以下の生成量のEAPで評価される
$ u_{a_2>a_1}^{(t)} \times u_{a_3>a_1}^{(t)} \times u_{a_1>a_4}^{(t)} \qquad (10.19)
$ 0.938となった
個々の命題が正しい確率の最小値以下の値となる
10.4. 推測例3(2水準の主効果の分析)
前節とは別人のS3氏が、自身で測定した「知覚時間」の実験のデータが表10-13
図10-3に平均値プロットを示し、表10-14に水準・交互作用の効果が$ 0より大きい、または小さい確率を示した
https://gyazo.com/a16f57e13833321aa3497dc903009baa
table: 表10-14 水準・交互作用の効果が0より大きい(小さい)確率 (S3)
a₁ a₂ a₃ a₄ b₁ (ab)₁₁ (ab)₂₁ (ab)₃₁ (ab)₄₁
0 < 0.564 0.310 0.387 0.731 0.997 0.292 0.489 0.650 0.575
≦ 0 0.436 0.690 0.613 0.269 0.003 0.708 0.511 0.350 0.425
$ b_1が$ 0より大きい確率は$ 0.997である
以上のことから、要因Bの主効果だけがあると判定される
水準が2つしかない要因の場合は、$ b_1が$ 0より大きい確率は、$ b_1が$ b_2より大きい確率と同義
したがって、$ b_1が$ b_2より大きい確率も$ 0.997
平均値プロットを観察すると「状態」の水準ごとの2本の折れ線が別の高さに描かれ、しかも2本とも大きな起伏がなく平坦
これが要因Bの主効果だけがある場合の特徴
表10-15に要因Bに相当する効果の大きさに関する指標を示す
table:表10-15 効果の大きさに関する生成量の推定結果 (S3)
EAP post.sd 2.5% 5% 50% 95% 97.5%
η²_b 0.090 0.054 0.008 0.015 0.084 0.190 0.212
δ_b 0.315 0.114 0.092 0.127 0.314 0.504 0.540
説明率は$ 0.090(0.054)[0.008, 0.212]
要因Bの主効果はデータの変動の$ 9.0\%を説明している
効果量は$ 0.315(0.114)[0.092, 0.540]
条件が一定に保たれているセル内での散らばりの$ 31.5\%の変動を要因Bの主効果が生み出している
放送授業
有意性検定の$ \alpha = 0.05には根拠がない
p値だけにたよって差の有無を判定することは誤り
差があるか否かという、実質科学的判定を、純粋に統計学の範囲内で済ませることはできない
差がある確率は、解釈可能だから、ドメイン知識を利用して確率を評価する